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We consider long simulations of 2D Kolmogorov turbulence body- forced by sinAyx on 
the torus {x,y) € [0,27r]^ with the purpose of extracting simple invariant sets or 'exact 
recurrent flows' embedded in this turbulence. Each recurrent flow represents a sustained 
closed cycle of dynamical processes which underpins the turbulence. These are used 
'"^ ■ to reconstruct the turbulence statistics in the spirit of Periodic Orbit Theory derived 

for certain types of low dimensional chaos. The approach is found to be reasonably 
successful at a low value of the forcing where the flow is close to but not fully in its 
asymptotic (strongly) turbulent regime. Here, a total of 50 recurrent flows are found 
with the majority buried in the part of phase space most populated by the turbulence 
giving rise to a good reproduction of the energy and dissipation probability density 
functions. However, at higher forcing amplitudes now in the asymptotic turbulent regime, 
I the generated turbulence data set proves insufficiently long to yield enough recurrent flows 

to make viable predictions. Despite this, the general approach seems promising providing 
enough simulation data is available since it is open to extensive automation and naturally 
generates dynamically important exact solutions for the flow. 

1. Introduction 

Ideas from dynamical systems have recently provided fresh insight into transitional 
and weak turbulent flows where the system size is smaller than the spatial correlation 
length. Viewing such flows as a trajectory through a phase space littered with invariant 
('exact') solutions and their stable and unstable manifolds has proved a fruitful way of 
understanding such flows (Eckhardt ct al. 2002, Kcrswell 2005, Eckhardt et al. 2007, 
Gibson et al. 2008, Cvitanovic & Gibson 2010, Kawahara et al. 2012). It is therefore 
natural to ask whether any ideas attempting to rationalise chaos may have something 
to say about developed turbulence. This is not to presuppose the two phenomena are 
simply related - that they are not has surely been appreciated for over 30 years - but 
merely an approach found useful in one may provide some insight into the other. One 
promising line of thinking in low-dimensional, hyperbolic dynamical systems stands out 
as a possibility - Periodic Orbit Theory. 

The study of periodic orbits as a tool to understand chaos has been a longstanding 
theme in dynamical systems dating back to Poincare's original work on the three body 
problem in the 1880s (Poincarc 1892; Ruclle 1978, Eckmann & Ruelle 1985, MacKay & 
Miess 1987) The fact that chaotic solutions can fleetingly, but also recurringly, resemble 
different periodic flows over time has always suggested that the statistics of the former 
may be expressible as a weighted sum of properties of the latter. However, this has 
generally remained a vague hope except for a special subclass of dynamical system where 
Periodic Orbit Theory has formalised this link (Auerbach et al. 1987, Cvitanovic 1988, 
Artuso et al. 1990a and for a recent review, Lan 2010). For these systems - very low 
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dimensional, fully hyperbolic invariant sets in which periodic orbits are dense ('axiom A' 
attractors) - there have been some notable successes (e.g. Artuso et al. 1990b, Cvitanovic 
1992 and later papers in the same journal issue, see also the evolving wcbbook Cvitanovic 
et al 2005). Here the invariant measure across the attractor can be expressed in terms of 
the periodic orbits which are dense within it so that crgodic averages can be determined 
from suitably weighted sums across the periodic orbits. Central to applying the approach 
is identifying a symbolic dynamics which can catalogue and order the infinite periodic 
orbits present in a chaotic attractor to give convergent expressions. 

Extending Periodic Orbit Theory to higher dimensional dynamical systems - most no- 
tably spatiotemporal systems - would obviously be highly desirable but represents a very 
considerable challenge. However, there are encouraging signs that something approaching 
this could be possible in fluid turbulence. The fact that a turbulent flow fleetingly yet 
recurringly resembles a series of smoother coherent structures or spatiotemporal patterns 
is a familiar observation perhaps first recorded by Leonardo da Vinci in his famous draw- 
ings and made mathematically by Hopf (Hopf 1948: see Robinson 1991, Holmes et al. 
1996 and Panton 1997 for overviews of subsequent work). Hopf 's vision of turbulence was 
of a flow exploring a repertoire of distinct spatiotemporal patterns where the implication 
was that these patterns were simple invariant solutions of the governing equations (e.g. 
equilibria, periodic orbits, tori etc - hereafter also referred to as recurrent flows). This 
viewpoint then advocates a dynamical systems approach even for fully turbulent flows. 
However, an attempt to build a prediction of turbulence statistics from the recurrent 
flows present is fraught with difficulties. Not only is there the daunting problem of ini- 
tially identifying enough of them in such high dimensional systems (typically 10^-10^ 
degrees of freedom) to make such a prediction seem feasible, but there is the problem of 
understanding how each should be weighted in any expansion. Finally, in the very likely 
eventuality that there is no symbolic dynamics for turbulence, it is impossible to know 
if important recurrent flows have been missed thereby compromising any prediction. 

The situation although very difficult, promises much and is not without hope. Efforts to 
extend the ideas of Periodic Orbit Theory to higher dimensional systems have focussed 
on 1-space and 1-time partial differential equations, most notably the 1-dimensional 
Kuramoto-Sivashinsky system (Christiansen et al. 1997, Zoldi & Greensidc 1998, Lan 
& Cvitanovic 2008, Cvitanovic et al. 2010) and the complex Ginzburg-Landau equation 
(Lopez et al. 2005). The emphasis in this work has mostly been to establish the feasibility 
of extracting recurrent flows directly from the 'turbulent' dynamics although some tenta- 
tive predictions were made (Christiansen et al. 1997, Lopez et al. 2005). The flrst attempt 
to extract a recurrent flow from 3 dimensional Navier-Stokes turbulence was made in a 
landmark calculation by Kawahara & Kida in 2001. In this work they managed to find 
one periodic orbit embedded in the turbulent attractor in a 15,422 degree-of- freedom 
(d.o.f.) simulation of small box plane Couette flow. This immediately raised the 'bar' of 
what had been thought possible and interestingly, they found that this one orbit was 
a very good proxy for their turbulence statistics. Van Veen et al. (2006) drew a similar 
conclusion albeit after discarding all but one of the few orbits they found when studying 
highly symmetric 3D body-forced box turbulence. Subsequent work in plane Couette flow 
by Viswanath (2007) essentially conflrmed the existence of Kawahara & Kida's (2001) 
periodic orbit (using 180,670 d.o.f.), found another and identifled 4 new relative periodic 
orbits (see also Lopez et al. 2005). These are periodic orbits where the flow repeats in 
time but drifts spatially in directions where the system has a continuous translational 
symmetry. Cvitanovic & Gibson (2010) report (using 61,506 d.o.f.) having identified 40 
periodic solutions, 15 relative periodic solutions with streamwise shifts and one relative 
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periodic orbit with a small spanwise shift in low Reynolds number and small box plane 
Couette flow. 

The state of the field is then that recurrent flows can be found in 3 dimensional 
Navier-Stokes turbulence calculations requiring up to O(fO^) d.o.f. (weak turbulence at 
low Reynolds numbers) but understanding how many can be found in a reasonable (toler- 
able) time and then identifying how dynamically important they are, remain outstanding 
issues. As a result, making useful predictions with any confidence using the set of recur- 
rent motions found seems some way off. With this background, our objective here is to 
make some contribution to this effort by mounting a systematic investigation of the issues 
in the simpler context of 2-dimensional Navier-Stokes turbulence. 

It's worth emphasizing that even if a 'turbulence' version of Periodic Orbit Theory 
ultimately proves beyond our grasp, the procedure of identifying different recurrent fiows 
buried within a turbulent solution has considerable value in its own right. This is because 
each recurrent fiow can be thought of as a sustainable dynamical process which helps 
underpin the turbulent state. Since they are 'closed' (recur exactly), their spatial and 
temporal structure can be dissected to reveal the fundamental physics involved. Just such 
an approach has helped uncover the 'self-sustaining process' ( Waleffe 1997) - stream- 
wise vortices generate streaks which are unstable to streamwise-dependent fiows which 
subsequently invigorate the streamwise vortices - in wall-bounded shear fiows following 
the discovery of a quasi-cycle in highly constrained plane Couette fiow by Hamilton, 
Kim & Waleffe (1995). Beautifully, this quasi-cycle turned out to indicate the presence 
of families of exact (unstable) travelling wave solutions to the Navier-Stokes equations 
(Waleffe 1997), the existence of which have revolutionized our thinking in transitional 
and weakly turbulent shear fiows (see the reviews by Kerswell 2005, Eckhardt et al. 2007 
and Kawahara et al. 2012). 

The specific framework investigated here is 2-dimensional Kolmogorov flow on a [0, 27r]^ 
torus (efficiently simulated using spectral methods) where the fiow is forced monochro- 
matically and steadily at a large length scale. This fiow has been extensively studied since 
Kolmogorov introduced the model in 1959 (Arnold & Meshalkin 1960) as a simple ex- 
ample of linear instability which could be studied analytically (Meshalkin & Sinai 1961). 
The fiow has many possible variations: torus aspect ratio (e.g. Marchioro 1986, Okamoto 
& Shoji 1993, Sarris et al. 2007), forcing wavelength (She 1988, Piatt et al. 1991, Arm- 
bruster et al. 1996), forcing form (e.g. Gotoh & Yamada 1986, Kim & Okamoto 2003), 
and 3-dimensionalisation (e.g. Borue & Orszag 1996, Shebalin & Woodruff 1997, Sarris et 
al. 2007). It has been experimentally realised using magnetohydrodynamic forcing (e.g. 
Bondarenko et al 1979, Obukhov 1983, Sommeria 1986) and latterly in soap films (e.g. 
Burgess et al. 1999). With an additional Coriolis term, Kolmogorov fiow can also be used 
as a barotropic ocean model on the /3-plane (e.g. Kazantsev 1998, 2001 and Tsang & 
Young 2008). 

The work by Kazantsev (1998, 2001) is particularly relevant for this study as this 
made the first attempt to apply Periodic Orbit Theory in a 211 d.o.f. discretization of 
a 2D Kolmogorov-like flow (differences include the addition of non-periodic boundary 
conditions, rotation and bottom friction). The work is most notable for his use of a min- 
imisation procedure to identify periodic orbits (59 found) as well as a good survey of 
relevant atmospheric literature. More recent work by Fazendeiro et al. (2010) (see also 
Boghosian et al. 2011) has started to study triply-periodic body-forced turbulence using 
Lattice-Boltzmann computations. Their focus was on developing another variational ap- 
proach for identifying periodic orbits based upon the idea of Lan & Cvitanovic (2004) 
and they describe convergence evidence for 2 periodic orbits. The approach starts with 
a closed orbit that does not satisfy the Navier-Stokes equations and uses a variational 
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method to adjust the orbit until it does. This requires manipulating the whole orbit 
at once and requires massive computations which are facilitated by the inherent paral- 
lelism of the Lattice-Boltzmann approach. In contrast, the approach adopted here is to 
start with a flow trajectory which does satisfy the Navier-Stokes equations but is not 
closed and to adjust the start of the the trajectory until it docs. This boils down to 
a Newton-Raphson root search in very high dimensions and iterative methods have to 
be employed to make things feasible. We adopt a Newton- GMRES-Hookstep procedure 
developed by Viswanath (2007,2009) and subsequently used with success by Cvitanovic 
& Gibson (2010) (see Duguet et al. 2008 for a slight variation which replaces the 'Hook 
step' with the 'Double Dogleg' step; Dennis & Schnabel 1996). 

The plan of the paper is as follows. Section 2 describes 2D Kolmogorov flow in detail, 
discusses its symmetries (§2.1) and makes connections with some previous direct numer- 
ical simulations (DNS) (§2.2). Key flow measures to be used subsequently are listed in 
§2.3. Section 3 describes the methodology used starting with the time-stepping code in 
§3.1, how initial guesses for recurrent flows are identified in §3.2, and then the Newton- 
GMRES-Hookstep algorithm in §3.3 (this draws its inspiration from Viswanath (2009)). 
§3.4 discusses how the algorithms were tested. Section 4 describes the results, first giving 
a flow orientation in §4.1, then reporting on how recurrent flows were actually extracted, 
before giving details of the recurrent flows found in §4.3. Section 5 describes an attempt to 
reproduce properties of 2D Kolmogorov turbulence before section 6 discusses the results 
and the outlook for future work. 



2. Formulation 

The incompressible Navier-Stokes equations with what is called 'Kolmogorov forcing' is 

— + It* • V*u* + - W = /^V*2tt* + xsin {2TTny*/Ly) x, (2.1) 

V*-tt*=0 (2.2) 

where p is the density, ly the kinematic viscosity, n an integer describing the scale of the 
(monochromatic) Kolmogorov forcing and x is the forcing amplitude per unit mass of 
fluid over a doubly-periodic domain [0, L^] x [0, Ly] (and in this section only * indicates a 
dimensional quantity). The system is non-dimensionalised by the lengthscale Ly/2'K and 
timescalc Ly/2-Kx so that the equations become 

du 1 

— +u-Vu + Vp= — V^M + sin(nj/) (2.3) 
ot Re 

V-u^O (2.4) 



where the Reynolds number is 



VX (Ly 



to be solved over the domain [0, 27r/a] x [0, 27r] (a := Ly/Lx). Given the doubly-periodic 
boundary conditions, dealing with the cross-plane vorticity equation is more natural and 
reduces simply to the scalar equation 

duj 1 

— — = z • Vx(mxwz) H V^cj — 71 cos (ny) (2-6) 

ot Re 

where loz := Vxm. (The form of the nonlinearity on the RHS is convenient for com- 
putation but can be further reduced to simply —u • 'Vu as the vortex stretching term 
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w. Vm = is, of course, absent in 2D.) Dealing with this equation is analogous to working 
with the streamfunction u = Vxi/'(.t, y)z since spatially-constant velocity and vorticity 
fields are not present so ip = V^^oj. 

2.1. Symmetries 

There is a shift-&-reflect symmetry 

7r 

S : [u,v]{x,y) -> [-u,v]{-x,y + -). (2.7) 

which shifts half a wavelength of the forcing function in y and reflects in x. Since there 
are n wavelengths in the domain, this transformation forms a cyclic group of order 2n — 1. 
There is also a rotation-through- tt symmetry 

'R:[u,v]{x,y)^[-u,^v]{-x,-y) (2.8) 

and the continous group of translations 

7i ■.[u,v\{x,y) ^ [u,v\{x ^l,y) for Q^l< — . (2.9) 

a 

The focus here is (unusually) not to take advantage of these, that is, the flow is allowed 
to fully explore phase space. 

2.2. Past literature 

Of all the previous work on 2D Kolmogorov flow, Piatt et al. (1991) seem to have carried 
out the most detailed study with n = A over the non-dimensional domain [0, 27r] x [0, 27r]. 
The same choices ?i = 4 and a = 1 were therefore made throughout the calculations 
reported here. With this. Re = and so the critical Reynolds number for linear 

instability is i?ec = 8v^ (i?ef'''" = %/2). Piatt et al. (1991) looked at the flow regime 
Re/Rcc ^ 3.54 over a 32 x 32 spatial grid so that 9.51 ^ Re ^ 33.6. Here we consider 
a 256 X 256 grid and look at 9.5 ^ iJe ^ 100 {Re/Rec ^ 10.5). Unfortunately, we were 
only able to confirm the detailed dynamics reported by Piatt et al. (1991) if we reduced 
our resolution down to theirs. 

The next closest study was She's (1988) which took = 8, a 64 x 64 grid and exam- 
ined y/2 Reshe 30 (26.9 < i?e sC 123.9 as Re = 8^/'^^/ Reshe) which corresponds 
to Re/Rec ^ 4.6. More recently, Sarris et al. (2007) considered 3D Kolmogorov flow 
over a variety of box aspect ratios considering 65 ^ ji^Sams ^ j^gQ jj^duding the 3- 
dimensionalisation of the flow considered here (then Re = SRe^'""^""). Typically they 
use 128 mesh points per wavelength of the forcing. At the time of writing, the world 
record for resolution when simulating doubly-periodic body-forced turbulence seems to 
be 32, 768^ (Boffetta k Musacchio 2010). 

2.3. Key measures of the flow 

Key measures of the flow which will aid the subsequent discussion are as follows (it = 
ux + vy): the mean flow, 

u{y) ■.= {u-x)t^,, (2.10) 

(initial conditions are such that {u{x,0) ■ y)^ = so that {u{x,t) ■ y)x — for all time 
); the bulk mean square of the fluctuations around the mean, 

^lnsit)-^{iu~uf)v, vlms{t)-^{^^)V: (2.11) 

and root mean square of the fluctuations as a function of y, 

Urmsiu) ■= J {{U - Uy)t,x, Vrmsiv) ■= J (w^)t,x ; (2-12) 
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the total kinetic energy and the kinetic energy of the fluctuation field 

E{t) := ^{u^)v, Et[t) := \{{u - (m)^)')^; (2-13) 
the total dissipation rate and the instantaneous power input 

:= ;^(|V«|2)v, I{t):={u&\n{ny))v, (2.14) 

with finally the laminar state, bulk laminar kinetic energy and bulk dissipation rate 
Re Re^ Re 

Ulam ■= ^ sin ny,X, Elam ■-^ -r^, Dlam'-=-^, (2-15) 

where the various averagings are defined as 

^ l'2-K / Q. 

( )x-= TT dx, 

27r Jo 

1 

{ )t lim - / dt. 



3. Methodology 

3.1. Time Stepping Code 

A 2D fully de-aliased pseudospectral code was used as developed in Bartello & Warn 
(1996). The original leapfrog+filter approach was replaced by the Crank- Nicolson method 
for the viscous terms and Heun's method (Euler predictor method) for the nonlinear and 
forcing terms so that only 1 state vector was required to accurately restart the code. 
This together with a constant time step size (except for the last step) means that the 
discretised flow is a dynamical system which closely matches the Navier-Stokes flow. 
Specifically, if ri(fc) = fft(u;(a;)) is the Fourier transform of w with k = {kx,ky), the 
vorticity equation (|2.6p in spectral space is 

^ = -Gn + f{n) (3.1) 

with G{k.x,ky):^^^^\^ (3.2) 
Ke He 

TL 

and f{n) := -i{ka;fft[uuj] + kyfft[vuj]) - -Sk^S(^\ky\-n)- (3.3) 

Here Si is the Kronecker delta function and takes the value 1 when i = and otherwise. 
A time step is performed by solving: 



At 2 

followed by solving: 



+ r!') + /(f!') (3.4) 
i(/(O^+i) + /(170) (3.5) 



At 2 

" V ' " ^ 

C-N Hcun 
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where the superscript is a time step index. With de-aliasing, a resolution of Nx x Ny 
corresponds in practice to the vorticity representation 

u:ix,y,t)= E ^Aty^"'''^'''^ (3.6) 

j = -Ar^/3/=-Af„/3 

where k = (aj,l), floo = and a mask is employed so that flji ~ for wavenumbers 
outside a specified domain S. Calculations reported here have a — I, n — A and = Ny 
so E := {{j, I) : j^ + P ^ is used. The number of active (real) degrees of freedom 

is therefore « which is « 22,800 (or exactly 22,428) for the = 256 used here 

(0 < j < 85 & -85 ; < 85 since n{-j, -I) = n*{j, I)). 

3.2. Near Recurrences 

The key idea pursued here is to extract recurrent flows directly from the turbulent DNS 
data with the implication that they are clearly dynamically important. With this in 
mind, the time stepping code was run for 10^ time units starting from random initial 
conditions and 'near recurrences' of the flow field searched for. These near recurrences 
were defined as episodes where 

uj{x + s,y + ^TTm,t + T) = uj{x,y,t) (S-''') 

'approximately' holds for some choice of the continuous shift ^ s < 2tt, the discrete shift 
m £ 0, 1, 2, . . . , n — 1 and T > over ^ x,y < 2tt. Periodic orbits correspond to s = 
m = and some period T > 0, travelling waves (TWs) to m = and s ~ cT with T > 
free where c is the phase speed, equilibria have s = m = and T free and relative periodic 
orbits have one or both of s and m not equal to zero with period T > 0. (The possible 
existence of relative periodic orbits - permitted by the inclusion of 2 free parameters (s 
and to) - is, of course, a reflection of the discrete and continuous translational symmetries 
of the system). The key is understanding what 'approximately' means, that is, how close 
should p.7p be to holding for it to signify the presence of a recurrent flow structure 
nearby. The only way to answer this seems to be to do computations and experiment. 
The search for near recurrences was done most efficiently by calculating every, say t=0.1 
or 0.2 steps, the normalised difference between states in wavenumber space suitably 
minimised over continuous shifts in x and discrete shifts in y as follows 

. . Ez |i^,7(0e-^-+^-wn _ ^^.^(^ , r)|2 

R{t, T) mm mm ir^ 3.8) 

0^.s<27rTne04,2,...,n-l EjEil%!(*)l 

where E, EJ%il' = a/{^nf j^'"''' ^^"^ uj'^dxdy, < Tthres ~ 0.5 < T < 100. Since 
R{t,0) = and dR{t,0)/dT > 0, the offset Tthres is defined adaptively as the first 
time at which dR{t,Tthres) / dT < 0. Figure [1] is a typical example of how R{t,T) looks 
as a function of t and T during a recurrent episode. The 9 black dots are the guesses 
identified by the code {R < Rthres = 0.3) over this time interval. All except one (the 
last dot at i w 171) subsequently converged to an exactly recurrent solution (the 4 dots 
for t < 130 to a periodic orbit (PI in Table 2) with period 5.3807 and the next 4 dots 
with t e [130, 160] to a TW (Tl in Table 2) with phase speed c = 0.0198). The threshold 
Rthres was choscn judiciously to give enough good quality guesses. 

3.3. UPO extraction method: Newton-GMRES-Hookstep 

Once a near-recurrence has been found by the above stated criterion, we then attempted 
to find if an exact recurrent flow was lurking nearby in phase space. This required a high- 
dimensional root finding algorithm acting on a state vector which completely specifies 




Figure 1: (Low resolution version for the arXiv) R{t,T) at Re = 40 contoured over 
t G [95,175] (x-axis) and T G [thres,50] (y-axis). The 9 black dots are the guesses 
identified by the code {R < Rthres) over this time interval. All except one (the last dot 
at t w 171) subsequently converged to exact solutions - the 4 dots for t < 130 to a 
periodic orbit with period 5.3807 (PI in Table 2) and the next 4 dots with t G [130, 160] 
to a TW (n in Table 2) with phase speed c = 0.0198 (i? values above 0.55 are not 
drawn/coloured for clarity. 



the velocity field 



X = 



s 
T 



(3.9) 



and contains information about the potential recurrence ($7 is a vector containing the 
scalars flji arranged in some fashion). The shift s is included since it can be adjusted 
continuously whereas the discrete shift m cannot and therefore is pre-set. To set up the 
Newton-Raphson algorithm (and we follow the excellent description in Viswanath 2009), 
it is convenient to define the infinitesimal generators and Ty of translations in a; and 



J{ajx+ly) 



2D Kolmogorov Turbulence 

y respectively 

Ny/3-l 

T^ui{x,y,t) ^ — = ^ ^ iajnji{t)e' 

-Ny/3 -N,/3 
Ny/3-1 N^/3-1 

-Ny/3 -N^/3 

as they act in spectral space 

Jl — > Jlx and Ty — > Jly 



did 

dx 

did 

dy 



(3.10) 



where each element ilji of O is mapped to iajilji in f2x and ilflji in Jly. Then, in spectral 
space the recurrence condition p.7p becomes 



F(r2o,s,T;m) := exp(.sT^ + i7rmTy)r2(f2o, T) - JIq = 



(3.11) 



where JIq = r2(t) and ft = n{t + T). If Xq = (Jig, so, To)"^ is an initial guess for a 
solution, then a better (next) guess Xq + Sxq = (fJo + S^l, sq + 5s, Tq + ST)'^ is given by 



Olio OS oT 



-F(f2o,so,7o;™) 



(3.12) 



These are dim(r2) equations for dim(Jl) + 2 unknowns. The extra two equations come 
from removing the degeneracy associated with these translational symmetries (the system 
is invariant under (x, t) — (x + s, t + T)). This can be done by imposing that Sfl, has no 
component which shifts the solution infinitesimally in the x-dircction or the direction 
(i.e. just redefines the time origin of the flow). The Newton- Raphson problem is then to 
solve 







dt 











sn 



5T 



F {no, So, To; m) 







(3.13) 



where Cls := exp(sTx + ^TTmTy)d is the 'back-shifted' final state and I is the dim(f2) x 
dim(J7) identity matrix. This is now in the standard form ASx = b with only the 
Jacobian matrix dCls/dCl not straightforward to evaluate {dCts/dT and dVlo/dt are 
found by substituting fig or rZg into the Navier-Stokes equations). 

Typically, the size of the matrix A is too large to store explicitly let alone attempt to 
solve A(5x = b directly. As a result, the only way to proceed is iteratively and GMRES 
(Saad & Schultz 1986) is convenient (see the excellent description by Trefethcn & Bau 
1997). Here only the effect of A on an arbitrary vector is needed. The effect of the 
troublesome Jacobian can be handled easily by a forward difference approach since 



OOs ^^s(»o + ey)-n,(no) 
dQ,o^ e 



(3.14) 
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where e is chosen such that ||ey|| = 10^^||r2o|| which balances truncation error with 
round-ofF error using double precision arithmetic and 1 1 • 1 1 is the Euclidean norm (using a 
more-physically orientated norm is clearly an interesting direction awaiting exploration). 

Straight Newton-GMRES is typically not good enough as guesses are usually not in 
the region where linearisation holds sufBciently well and divergence to infinity is com- 
monplace. Instead it proves useful to modify the approach to incorporate a trust region. 
Following Viswanath (2007, 2009), we use the 'hook-step' method (Dennis & Schnabel 
1996 §6.4.1) which can be easily built on top of the GMRES process. Exactly how the 
approach is implemented can vary and we adopt what looks to be slightly different algo- 
rithm to Viswanath (2007, 2009) in which GMRES is used first to derive an approximate 
solution to A^x = b before this 'solution' dx is moved into the trust region. The ad- 
vantage of this is there is a clear convergence criterion that can be imposed to terminate 
the initial GMRES algorithm. Before stating this, it's worth first briefiy describing the 
GMRES algorithm itself which is based upon a simple idea. The GMRES algorithm for 
solving A^x = b at iteration n approximates <5x by the vector <5x„ in the Krylov space 
3Cn '■= ( b, Ab, A^b, . . . , A"^-'-b ) that minimises the norm of the residual 

||A^x„ - b|| (3.15) 

(Trefethen & Bau 1997). For numerical stability, an orthonormal basis for is con- 
structed using a Gram-Schmidt-style iteration as follows 

9i = > = Aq,, Qi+i = q^+l/\\q^+l\\ iel...,n-l 

(3.16) 

so that if Q„ is the matrix with columns Qi, 52, • • • , g„, then AQ„ = Q„+iH„_|_i_„ where 
H„_|_i^„ is the upper (n + l) x n left section of an upper Hessenberg matrix generated by 
the basis orthonormalisation (Trefethen & Bau 1997, p252). With this basis the solution 
Sku = QnYn (where y„ is an n- vector) minimises ||Qn-i-iH„+i,„y„ — b|j {N equations 
and n unknowns) or equivalently ||H„+i.„y — ||b||ei|| (n+l equations and n unknowns) 
since the only non-zero entry in Q^_^j^b is the first entry. This can be accomplished by 
a singular value decomposition (SVD) of H„+i^„ into U„+iDV^ (where U„+i and V„ 
are orthonormal matrices and D is an (n -I- 1) x ri diagonal matrix with a zeroed bottom 
row) through straightforwardly solving the first n equations Dz„ = p„+i := ||b||U^_|_]^ei 
followed by V^y„ = z„ and then ^x„ = QnYn - The modulus of the remaining unbalanced 
component pn+i{n + 1) then gives the minimum value or residual. The iterations are 
continued until 

||A^x„-b|| _ ||Dz„-p„+i|| _ |p„+i(n+l)| ^^^^ ^^^^^ 



\m \\Pn+l\\ \\Pn+l\ 

where tol is a small number typically chosen in the range 10""* to 10~^ (the majority of 
the computations reported here were obtained using a value of 10"'^). If ||F(x -f- (5x„)|| 
is not smaller than ||r(x)|j or more specifically not well predicted by the linearisation 
around x i.e. ||F(x) + A<5x„||, then the approximate solution of the linearised problem 
is transformed back to a smaller trust region where the linearised problem is valid. This 
is done by adding the constraint ||<5x„|| ^ A or equivalently ||y„|| ^ A to the GMRES 
minimisation p.lSp : this is the hook step. The beauty of this adjustment is that it is 
a very natural modification of the GMRES approximate solution since the (innermost) 
problem for is then 



min||Dz„ - p„+i|| s.t. ||y„j| = ||zn|| ^ A 



(3.18) 
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Constructing the Lagrangian 

L := (Dz„ - p„+i)2 + f^izl + A^) (3.19) 

where /z is a Lagrange multiplier imposing the trust region constraint, leads to the min- 
imisation equations 

2d{i){d{i)zrS) - Pn+i(«)) + 2nznii) = (3.20) 

2/i/? = (3.21) 
zl + - = (3.22) 

where di is the ith diagonal clement of D. The solution to this is 

,^^^) :^ ^^flM l^^^n (3.23) 
O/^ I /z 

with either = as ||z„|l < A (the original GMRES solution) or yu 7^ chosen so that 
||z„|| = A (in practice fi is just increased until ||z„|| < A). An acceptable solution (5x„ 
is signalled by 

F (x + ^x„ ) F (x) + cA^x„ (3.24) 

where some value of c € (0, 0.5) is chosen (eqn 6.4.14, Dennis & Schnabel 1996: we took 
the least demanding value of c = 0). If this does not hold, A is decreased and the hook 
step repeated until it is. Depending on how easily this improvement condition is met, the 
trust region may be relaxed (e.g. if linearisation holds well - F(x + Jx„) sa F(x) + AJx„) 
or not for subsequent Newton steps. 

This algorithm can be readily extended to perform solution branch continuation: see 
appendix A. Furthermore, since we know how to calculate the action of the Jacobian on 
any vector (sec (|3.14p . the linear stability of an exactly recurrent flow can also be readily 
found using the Arnoldi technique (e.g. using ARPACK to extract extremal eigenvalues). 

3.4. Testing 

The modified (Crank-Nicholson+Hcun) time-stepping code was thoroughly validated 
against the well-tested Leapfrog+filtcr code developed by BartcUo & Warn (1996). The 
Newton-GMRES-Hookstep algorithm developed on top of this was tested by attempting 
to converge onto a known periodic orbit. This orbit was originally found by tracing bifur- 
cations up from the basic state. For n = 4, the ID basic state becomes linearly unstable 
at Re = 9.9669 for disturbances 27r-periodic in x giving rise to a steady 2D state which is 
3l-symmetric. This state loses stability to a stable periodic orbit within the 3l-symmctric 
subspace for 30 < Re < 31 before this orbit becomes unstable at Re > 32 through a 
torus bifurcation. The periodic orbit at Re = 31 was easily found by time-stepping within 
the 3^-symmctric subspace yet is unstable to 3^-asymmctric disturbances in the full un- 
restricted space. Having such an orbit to experiment with was invaluable for building up 
confidence in the code and some feel for how the tolerances of the algorithm should be 
set (e.g. tol in ((3Tf|) V 

4. Results 

4.1. Flow Orientation 

2D Kolmogorov flow is linearly unstable at a comparatively low Re which depends 
strongly on the imposed periodicity in the forcing direction: see figure [21 For the domain 
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Figure 2: The energy stability Reynolds number Res (red dashed line) and the linear 
instability Reynolds number Reun (blue solid line) for sin4y£ forcing over the torus 
[0,27r/a] X [0, 27r]. Note that Reun — oo as a — ?> 4 so the domain is squeezed down 
to [0,7r/2] X [0, 27r]. The dotted line at a = 1 is the present case {Rce = 6.8297 and 
Reiin = 9.9669). Res ^ and Reun ^ « 9.5137 as a ^ 0. 



studied here [a — 1), disturbances to the base flow (|2.15|) fail to decay monotonically at 
ReE = 6.8297 and then start to grow exponentially at Reun ~ 9.9669. Figure [3] shows 
that this initial bifurcation is to a steady flow {D/Diam < 1 and Et/E = 0) until Re « 15 
whereupon time dependence appears. For 15 < Re < 23, some metastability is noticed 
which is illustrated in figure|3]at Re = 22 for two different initial conditions. One leads to 
a chaotic- looking dissipation signal across the time interval [500, 1500] whereas the other 
drops out of this chaotic state at just over t=1000 to converge on a stable travelling 
wave solution (later named Tl). Beyond Re ~ 23, the chaotic state presumably becomes 
an attractor or the probability of dropping out of this state becomes so small that it is 
not picked up over the time windows studied (10^ units here and 10^ later). Finally, an 
asymptotic regime is approached for Re > 50. The preliminary calculations performed 
here for 100 < Re < 200 tentatively support the asymptotic scaling laws D ^ i?e~^/^ 
and U \/2E - Re'^/^ although the noisy data clearly warrants much longer time 
averaging to confirm this. 

Given this general flow behaviour, we chose to concentrate on analysing the flow at 
Re = 40 (approaching the asymptotic regime), and three values. Re — 60,80 and 100, 
which get deeper into the asymptotic regime. Figures [Hand Ogive an idea of the temporal 
and spatial scales in the flow at the two extremes. Re — 40 and Re = 100, of our 
study. Both indicate a hierarchy of temporal and spatial scales (which broaden with 
Re) indicative of 2D turbulence. Figure [6] confirms that the flows studied for Re ^ 100 
are well-resolved: there is 10 orders of drop off in the enstrophy spectrum in the most 
demanding case {Re = 100) for the 256^^ resolution used throughout this work. 
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200 



Figure 3: The normalised dissipation D/Diam (blue solid line with dots) and fractional 
fluctuation kinetic energy Et/E (red dotted line connecting dots) to illustrate how the 
flow changes with Re. The first bifurcation at Reiin = 9.9669 is steady so D/Diam 
decreases below 1 but Et/E remains until time-dependent fiows appear at around 
Re = 15. Curves were generated by using a random initial condition, running for 500 
time units and then calculating averages over the next 1000 time units. This produces a 
unique average except over the interval 15 < Re < 23 where there is some metastability 
as shown. Inset: Re^^^D (upper blue line) and iRe^^/^E^^^ (lower thick red line) against 
Re on a log scale showing the apparent start of the asymptotic regime D ^ Re~^^^ and 
U := \/2E ^ Re^/^ . The downward sloping lowest line shows U jRe^l"^ for comparison. 



4.2. Finding recurrent structures 

A standard hunt for recurrent fiows involved integrating the flow from random initial 
conditions for a period of 10^ time units. Initially, Rthres was set at 0.15 for three runs 
at Re = 40 (labelled a, 6 and c in Table 1) which produced only 9, 7 and 13 guesses 
respectively. Relaxing Rthres to 0.3 (run d), however, produced 885. This threshold value 
proved adequate at Re = 60 (nearly 300 near-recurrences detected over runs e,/ and g) 
but had to be further relaxed to 0.35 at Re = 80 and 0.4 for Re = 100: see Table 1. 
Unfortunately, it was noticed after these (Series A) runs had been completed and the 
guesses tested for convergence that only s = m = shifts had been searched over. So 
the runs were repeated (Series B runs o, p, q and r) searching specifically for recurrences 
which selected either s 7^ and/or m 7^ to minimise R. This was done to indicate 
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1500 



Figure 4: The normalised dissipation D{t)/Diam vs t € [500,1500] for Re=22 (two dif- 
ferent initial conditions: top solid blue & green dashed lines), Re = 40 (second lowest 
red line) and Re = 100 (lowest black line). Re = 22 shows chaotic saddle behaviour with 
one trajectory seemingly dropping out back randomly to the travelling wave Tl (stable 
at Re = 22). 



the frequency of observing strictly periodic near-recurrences and relative periodic near- 
recurrences. 

The initial trawl for near-recurrences took a few weeks (each case run on a Xeon 
X5670 processor) with the DNS code slowed considerably by the need to search for near- 
recurrences every 0.1 or 0.2 units in time (which is anything from 20 to 100 numerical 
time steps). The more time-consuming activity, however, was attempting to converge 
the near-recurrent guesses to exact solutions. Adopting fairly conservative limits for the 
Newton-GMRES-Hookstep procedure - maximum period considered was 100, maximum 
number of Newton, GMRES and Hook steps were 75, 500 and 50 respectively - typically 
lead to run times of a couple of months for each of the Re = 60,80 and 100 runs. The 
data for Re = 40 had to be subdivided 12 ways to make the process manageable. These 
numbers make it clear why a very efficient DNS code was important for this work. 

Table 1 also indicates the conversion rate of near-recurrences guesses to exactly recur- 
rent solutions. There is considerable duplication of such solutions so that a much smaller 
set of distinct recurrent structures is obtained. 
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Figure 5: A time sequence of vorticity plots at Re ~ 40 (upper) and Re = 100 (lower) 
over the 2t: x 27r domain. Time difference between plots is 1 time unit and 10 contours 
arc drawn between the maximum (white/light) and minimum values (red/dark) of the 
perturbation vorticity -10.25 ^ w s$ 8.6 at Re = 40 and -18.4 s$ cj 20.3 at Re = 100. 



Figure 6: The enstrophy spectrum for the rightmost snapshots in figure[5]at Re ~ 40 (blue 
circles) and Re = 100 (red crosses) D(n) := V i — - i |ri,;P for n — 1,2. ..85: 

II— 2^!vj +' <"+2 

recall definition p. 61 with a — 1) (e.g. there are 264 and 516 wavenumbcrs included for 
n = 40 and 84 respectively). The energy spectrum E{n) w D{n)/n^ has a steeper drop 
off. The dotted line indicates the wavenumber of energy injection. 



16 Gary J. Chandler and Rich R. Kerswell 





Run 




dt 


duration 


^ of guesses 


^ of convergences 


Series A 


(s = m = 0) 












Re = 40 


a 


0.15 


0.005 


10' 


9 


5 




b 


0.15 


0.005 


lO' 


7 


3 




c 


0.15 


0.005 


10' 


13 


5 




d 


0.30 


0.005 


10' 


885 


553 


Re = 60 


e 


0.30 


0.003 


10' 


102 


64 




f 


0.30 


0.003 


10' 


104 


67 




g 


0.30 


0.003 


10' 


78 


58 


Re = 80 


h 


0.35 


0.0025 


10' 


53 


31 




i 


0.35 


0.0025 


10' 


60 


37 




j 


0.35 


0.0025 


10' 


41 


25 


ite = iOO 


1 


0.4 


0.002 


iO 


to 


Q A 

34 




m 


0.4 


0.002 


10' 


91 


33 




n 


0.4 


0.002 


10' 


93 


42 


Series B 


(s / or m / 0) 












Re = 40 


o 


0.30 


0.005 


10' 


1223 


540 


Re = 60 


P 


0.30 


0.003 


3 X 10' 


163 


7 


i?e = 80 


q 


0.35 


0.0025 


3 X 10' 


66 


15 


Re = 100 


r 


0.4 


0.002 


3 X 10' 


84 


12 



Table 1: DNS data used to extract equilibria, travelling waves, periodic orbits and rel- 
ative periodic orbits at Re = 40,60,80 and 100. The value of Rthres cannot be set too 
ambitiously. Run d yielded all the solutions thrown up by runs a,b and c combined. 



4.3. Recurrent structures found 

4.3.1. Re=40 

Table 2 lists the recurrent structures found at Re = 40. The equilibrium flow El 
(see figure [7]), which was found many times in the series A runs, corresponds to the Jl- 
symmetric steady state which bifurcates off the basic solution at Re = 9.9669 as shown 
in figure [H El loses stability at about Re = 15 to the travelling wave Tl or later via PI 
in the 3?-symmetric subspace for a Re G (30, 31) {P2 which is !K-symmetric and P3 which 
is not bifurcate at yet higher Re from El). All three flows. El, PI and Tl are found 
to be repeatedly visited by the (series A) DNS indicating the strong influence of the 
3l-symmctric subspace on the 'turbulent' dynamics despite them all being unstable (e.g. 
El has 9 unstable directions at Re = 40; see Table 2). However, a further 47 recurrent 
flows were also identified from the DNS: another travelling wave T2, two further periodic 
orbits P2 and P3, and 44 relative periodic orbits, Rl — R6 and R18 — R55 (note all 
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UPO frequency 



T 



N 



Er=i»e(A,)(max»e(A,)) 



El 


> 261 











9 


1.296(0.249) 


Tl 


127 


0.0198 









4 


0.142(0.068) 


T2 


1 


0.0096 









4 


1.227(0.454) 


PI 


143 




5.380 







7 


0.570(0.191) 


P2 


6 




2.830 







5 


0.742(0.223) 


P3 


2 




2.917 







7 


0.992(0.236) 


Rl 


1 




56.677 


0.092 





3 


0.156(0.077) 


R2 


1 




25.401 


0.199 





5 


0.254(0.123) 


R3 


1 




54.280 


0.200 





3 


0.195(0.108) 


R4 


1 




6.720 


0.106 





8 


0.870(0.343) 


R5 


1 




23.780 


0.022 





4 


0.376(0.156) 


R6 


4 




20.808 


0.060 





3 


0.258(0.172) 


R18 


1 




37.233 


0.270 





5 


0.242(0.165) 


R19 






12.207 


0.243 





2 


0.141(0.070) 


R20 






16.586 


5.827 


1 


4 


0.289(0.103) 


R21 






17.470 


5.765 


3 


5 


0.348(0.143) 


R22 






19.723 


0.222 





4 


0.297(0.172) 


R23 






19.762 


0.513 





4 


0.302(0.127) 


R24 






19.779 


6.035 





5 


0.292(0.202) 


R25 






20.201 


5.898 


3 


6 


0.380(0.138) 


R26 






20.385 


1.334 


2 


7 


0.714(0.270) 


R27 






20.632 


5.871 


3 


4 


0.365(0.127) 


R28 






20.885 


5.987 


1 


4 


0.360(0.121) 


R29 






20.909 


0.306 


1 


5 


0.380(0.124) 


R30 






21.310 


5.694 





5 


0.330(0.100) 


R31 






21.725 


5.799 





3 


0.319(0.133) 


R32 






22.560 


0.006 


1 


4 


0.283(0.096) 


R33 






22.617 


5.660 





5 


0.478(0.156) 


R34 






23.157 


0.265 





3 


0.260(0.113) 


R35 






23.417 


5.936 


3 


4 


0.489(0.183) 


R36 






24.465 


6.010 


3 


4 


0.358(0.191) 


R37 






25.870 


0.182 





3 


0.272(0.122) 


R38 






25.934 


0.227 





4 


0.263(0.125) 


R39 






27.138 


6.248 





4 


0.391(0.107) 


R40 






28.817 


5.971 





5 


0.238(0.116) 


R41 






32.541 


0.349 


1 


4 


0.224(0.153) 


R42 






34.316 


5.886 





5 


0.163(0.120) 


R43 






34.530 


5.742 





3 


0.220(0.134) 


R44 






34.917 


-0.059 


3 


4 


0.325(0.139) 


R45 






36.549 


6.027 


3 


3 


0.183(0.118) 


R46 






36.627 


0.197 


3 


4 


0.202(0.139) 


R47 






36.812 


5.648 





3 


0.155(0.074) 


R48 






37.079 


6.103 





4 


0.171(0.134) 


R49 






37.233 


0.270 





3 


0.241(0.165) 


R50 






37.698 


3.499 


1 


6 


0.477(0.146) 


R51 






39.368 


6.070 





5 


0.192(0.098) 


R52 






39.619 


0.380 





5 


0.242(0.067) 


R53 






41.400 


5.806 





4 


0.176(0.081) 


R54 






49.645 


6.054 


3 


4 


0.160(0.065) 


R55 






53.073 


6.031 





4 


0.189(0.105) 



Table 2: All the invariant sets found directly from turbulent DNS data (from series A 
above the separating blank line and from series B below). 'Frequency' is the number of 
times the solution was extracted (Series A runs). There is one steady Equilibrium, c is 
the phase speed of the Travelling waves found, T is the Period of periodic and Relative 
periodic orbits which also either have a shift s and/or shift m. N is the number of unstable 
directions and Xl^i ^^i^j) the sum of the real parts of all the unstable eigenvalues. 
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UPO frequency c T -s -m N J^;^'!^ 3ffe(A^) ( max 5Re(Aj) ) 



El 


4 











14 


5.053(0.858) 


Tl 


high 


0.0019 









4 


0.139(0.064) 


T3 


high 


0.0124 









17 


3.377(0.684) 


T4 


1 


0.0082 









3 


0.257(0.178) 


R7 


high 




2.472 


0.036 





9 


0.911(0.214) 


R8 


1 




1.638 


0.022 





14 


2.903(0.681) 


R56 






16.326 


0.588 


2 


6 


0.609(0.139) 


R57 






17.909 


5.802 





7 


0.805(0.169) 


R58 






20.546 


0.659 


2 


8 


0.529(0.168) 


Tl 


high 


0.0115 









6 


0.360(0.105) 


T3 


8 


0.0154 









21 


5.588(0.958) 


T5 


3 


0.0831 









20 


4.183(0.658) 


R7 


10 




2.299 


0.054 





13 


1.326(0.181) 


R8 


2 




1.705 


0.028 





18 


4.318(1.105) 


R9 


1 




2.150 


0.032 





19 


3.987(1.026) 


RIO 


1 




1.280 


0.020 





20 


5.310(0.878) 


Rll 


1 




2.443 


0.031 





10 


0.704(0.277) 


R12 


1 




2.095 


0.034 





11 


3.083(1.004) 


R13 


1 




15.285 


0.181 





8 


0.409(0.131) 


R59 






15.667 


0.397 


1 


11 


0.949(0.176) 


R60 






16.071 


0.462 


1 


11 


1.136(0.248) 


Tl 


high 


0.0155 









10 


0.646(0.122) 


T3 


7 


0.0179 









25 


6.491(1.042) 


T4 


5 


0.0118 









3 


0.684(0.370) 


T5 


6 


0.0691 









28 


6.238(0.689) 


P4 


1 




1.185 







16 


7.376(1.201) 


R7 


4 




1.971 


0.030 





15 


1.933(0.326) 


Rll 


3 




2.262 


0.001 





9 


0.905(0.385) 


R12 


1 




1.902 


0.029 





14 


3.953(1.244) 


R14 


5 




4.526 


0.071 





8 


0.428(0.105) 


R15 


2 




1.984 


0.122 





16 


3.110(0.556) 


R16 


2 




1.938 


0.121 





6 


0.945(0.270) 


R17 


1 




3.827 


0.008 





16 


2.894(0.818) 


R61 


1 




1.344 


0.090 





21 


4.997(0.476) 



Table 3: All the invariant sets found directly from turbulent DNS data (from series A 
above the separating blank line and from series B below: none were found at Re = 100 
in the series B runs despite the numerical gap between R17 and i?61). 'Frequency' is 
the number of times the solution was extracted. Solutions listed under each Re indicate 
those actually extracted at that Re - hence multiple entries. 'Frequency' is the number 
of times the solution was extracted (Series A runs). There is one steady Equilibrium, 
c is the phase speed of the Travelling waves found, T is the Period of periodic and 
Relative periodic orbits which also either have a shift s and/or shift m. N is the number 
of unstable directions and J2f=i^^i^j) ^^"1 of the real parts of all the unstable 
eigenvalues. 
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have a non-zero shift s and some also a non-zero integer m). A priori, we expected to 
find mainly small period recurrent structures due to the method of extraction. Longer 
periods mean more time for the turbulent trajectory to diverge away from the unstable 
recurrent flow and hence a higher probability for a) the episode to escape detection as a 
nearly recurrent flow and b) even if detected, for GMRES to fail to converge due to the 
quality of the initial approximation. This seems borne out by the periodic orbits found 
but not for the relative period orbits where the majority have a period over 20 and some 
over 50 time units. That such long period structures exist and were 'extractable' from 
the DNS frankly was a surprise and begs the question whether our 'long' runs of 10^ 
time units (now known to be only a factor of 0(1000) longer than some recurrent flows) 
were actually really long enough to capture all the structures possible. This issue will be 
raised again below. 

With so many recurrent flows found, it becomes impractical to display and characterise 
each flow separately. Table 2 lists some key characteristics along with their stability in- 
formation (all are unstable but none with more than 9 unstable directions out of 22,428 
possible directions). One useful projection used by Kawahara & Kida (2001), however, 
is the 'energy out {D{t)) verses energy in (/(i))' plot which is shown in figure [5] (both 
quantities normalised by Diam). The line D — I corresponds to dissipation exactly bal- 
ancing energy input which has to be the case over all times for equilibria and travelling 
waves (which are just equilibria in an appropriate Galilean frame): these are therefore 
just points on this line in this plot. Figure [8] shows how a representative subset of these 
recurrent fiows look when compared with the joint dissipation-input probability den- 
sity function (pdf) of the DNS. The darkest shading makes it clear that the DNS stays 
predominantly in the region 0.055 < //Aam < 0.115,0.06 < D/Diam < 0.11. The re- 
current fiows shown are also dominantly concentrated in this region although there are 
two relative periodic orbits shown - i?26 and i?50 - which have large dissipation episodes 
(it's worth emphasizing that the basic state would be represented by the point (1,1) in 
this plot so the turbulent flow adopts a much reduced dissipative state). Since this D 
verses / plot is such a drastic projection of the dynamics, the fact that two flows look 
close there doesn't necessarily mean they are close in the full phase space. However, 
because all the recurrent flows discussed here have been extracted from turbulent DNS, 
this conclusion nevertheless seems reasonable. 

In figure |9] we focus on one typical 'embedded' relative periodic orbit R25 which stays 
within the central region of the DNS joint pdf. There is a clear temporal cycle where the 
energy input increases (exceeding the dissipation) and then decreases (now exceeded by 
the dissipation). Plotting the associated vorticity fields over this cycle - figure [10]- shows 
the character of the fiow. At the dissipation low point (time 17 in figures l9l and ITOl) . the 
vorticity is concentrated into weak y-aligned patches which are separated from each other 
whereas at the high dissipation point (time 8), the vorticity seems to be undergoing a 
shearing episode with only one stronger vortex recognisable. These two extremes bear 
more than a passing resemble to either Tl [D/Diam = 0.071) or T2 (0.071) and El 
(0.102) respectively suggesting that i?25 is probably a closed trajectory linking their 
neighbourhoods. i?50, in contrast, undergoes a large high dissipation excursion as shown 
more completely in figure [TT] The associated vorticity fields - see figure [12]- show similar 
structures to i?25 when in the same part of (/, D) space (compare t=5 for i?25 with t=0 
for i?50, and t=15 for i?25 and t=31 for i?50) but i?50 exhibits intense shearing too and 
vortex break-up at times 5, 8 and 9. i?50 clearly reflects an important but infrequent 
aspect of the turbulent dynamics as indicated by the fact that the joint (D, /) pdf of the 
DNS stretches to such high values of the dissipation. Whether we have extracted enough 
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Figure 7: The steady solution El (left) and the travelling waves Tl (middle) and T2 
(right) at Re = 40. Vorticity is contoured using 15 contours between -7.2 (dark, red) 
and 12.1 (light,white) (range is -6.5 ^ w < 12.1 for El, -7.2 ^ w < 7.2 for Tl and 
-6.8 w =^ 10.4 for T2). 

of such recurrent structures to capture this episodic behaviour is of course a key issue 
for this study and will be discussed in 

Figure [T2] is an attempt to show more of the recurrent structures found by zooming 
in on the central dashed box drawn in figure [51 This illustrates the intricacy of most of 
the flows found - many of the relative periodic orbits trace complicated D ~ I curves 
whereas, in contrast, the periodic orbits are simple loops. Another key observation is that 
some relative periodic orbits look very similar - e.g. R28 and R29 (and other pairings not 
shown). This, of course, resonates with the mental picture one has of periodic orbits being 
dense in a chaotic attractor. In fact, the consecutive numbering of R28 and R29 indicates 
that these relative periodic orbits were found concurrently from the DNS confirming their 
proximity in phase space. Also it is clear than some relative periodic orbits look like 
merged versions of two shorter orbits (not shown) again consistent with low-dimensional 
dynamical systems thinking. 

4.3.2. Re^60, 80 & 100 

At higher i?e, we managed to extract far fewer recurrent flows from the DNS. There 
are certainly reasons to expect this, most notably that the recurrent flows present should 
become more unstable and it is therefore harder to find good guesses from the DNS. 
And there is also the fact that the 'turbulence' should explore more of phase space and 
therefore close visits to simple invariant sets should become rarer. However, the sharp 
drop in the number of recurrent flows found - see Table 3 - was still a surprise. In keeping 
with the philosophy of this work, only recurrent flows extracted from the DNS at that 
Re are listed in Table 3 under the relevant Re heading. This then says nothing about 
whether a certain recurrent flow found at one Re might not exist at another. To explore 
this a little, we carried out some branch continuation (see appendix A) on the recurrent 
flows extracted from the series A DNS while the runs and analysis for the series B DNS 
were progressing. The results are shown in figure [14] colour coded to group recurrent 
flows found at the same Re and with black dots indicating branches detected at a given 
Re (note the rescaled dissipation measure on the ordinate to make the plot clearer). For 
example, the r3 branch is shown as a dashed red line (second red dashed line up from 
the Re axis) as it was first found at Re = 60 and it bears three black dots marked at 
Re = 60, 80 and 100 as T3 was extracted from the DNS at all three Re. This type of 
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Figure 8: The normalised dissipation D{t) / Diam verses I / Diam for a small collection of 
the recmTcnt flows found with the pdf of the DNS turbulence plotted in the background 
(11 shades at levels 10" where a ^ -5, -4.5, . . . , -0.5, 0). Plotted are El {D/Diam = 
I/Diam = 0.102), Tl (0.071), T2 (0.070), PI, P2, Rl - 6, R18, R26 (large blue orbit) 
and i?50 (the even larger magenta orbit). The dashed box is used to show more recurrent 
structures in figure [T3l 

analysis can indicate the bifurcation structure - e.g. T3 clearly bifurcates off a recurrent 
flow found at Re = 40 - but is very time-consuming to pursue through to completion as 
branches can become difficult to continue and interpret (note the number of open circles 
in figure [14] which indicate where the branch continuation procedure stagnated for some 
reason). This aside, the overriding impression is one of simple invariant sets proliferating 
with increasing Re. Notably, only two recurrent flows found at Re = 40 are also extracted 
from the Re = 60 DNS - El (the highest blue line with a dot at Re = 60) and Tl. El 
seems to lose dyanmieal importance for yet higher Re but Tl is found for all four Re 
studied here. 

Figure [15] shows the new travelling waves found and figures [16] to [18] the D — I plots 
where now all the recurrent flows found at the respective Re are marked. Again most 
sit in the D — I region where the DNS spends the majority of its time although as at 
Re = 40 there are some outliers (e.g. El at Re = 60, R8 at Re = 80, and PA and 
RIA at Re = 100). That i?14 actually appears outside the footprint of the DNS pdf 
at first looks erroneous but is in fact merely an indication that when the 'turbulence' 
approached i?14 in phase space, it maintained higher (global) dissipation and energy 
input than i?14. This can occur when part of the domain resembles i?14 while the rest 
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Figure 9: The normalised dissipation D{t) / Diam verses I / Diam for the relative periodic 
orbit R25 at Re=:40 with the pdl of the DNS turbulence plotted in the background (11 
shades at levels 10" where a — —5, —4.5, . . . , —0.5, 0). The labels refer to times along the 
orbit at which snapshots are shown in figure [TU) 



does not and exhibits enhanced dissipation. A good example of this is shown in figure 
[19] which details the turbulent episode which signalled the presence of P4 (left column) 
alongside the successfully converged periodic orbit P4 (right column). Visually, the eye 
is drawn to the centre of the domain where in both columns an isolated vortex is clearly 
seen rotating in a clockwise fashion. However, the corners are just as significant in that 
they also indicate an isolated vortex, yet this is stronger with higher gradients (and 
hence larger dissipation) for P4 than the DNS signal. Plotting the two time sequences 
on a D — I plot shows P4 as a closed loop much higher up the D = I line than the DNS 
(not shown). 



5. Recurrent flows as a turbulent alphabet 

Given the sets of recurrent fiows extracted at each Re, the question is then how to 
use them to predict properties of the turbulence encountered. Periodic Orbit Theory 
advocates a weighted expansion of the recurrent flows such that 



prediction 



T„=i Wi 



(5.1) 
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Figure 10: A time sequenee of vortieity plots for i?25 at Re = 40 at times (running left 
to right across the top and then bottom) t = 0,5, 8, 10, 13, 15, 17, 20 (marked as dots on 
figure [9]). The period is 20.2 so the flow in the bottom right is nearly the same as the top 
left except for shifts in x and y. For all, 15 contours are plotted from -12 to 12. 
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Figure 11: The normahsed dissipation D{t)/Diam verses I/Diam for the relative periodic 
orbit i?50 at Rc=40 with the pdf of the DNS turbulence plotted in the background (11 
shades at levels 10" where a = —5, —4.5, . . . , —0.5, 0). The labels refer to times along the 
orbit at which snapshots are shown in figure 1121 
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Figure 12: A time sequence of vorticity plots for i?50 at Re = 40 at times (running left 
to right across the top and then bottom) t = 0,5,8,9,10,26,31,37 (marked as dots on 
figure [TT|l . The period is 37.7 so the flow in the bottom right is nearly the same as the 
top left except for shifts in x and y. For all, 15 contours are plotted from -12 to 12. 



where F is any property such the mean dissipation rate, the mean profile or a pdf, N the 
total number of recurrent flows and the weights are 



Wi 



n 1(1 



e 



(5.2) 



(e.g. Cvitanovic 1995, Lan 2010). Here A^,'-* is the fcth eigenvalue of the linearised operator 
around the ith recurrent flow of period T*^'^ (the associated Floquet multiplier would be 

l^k^ := C^fc'^'"'). Calculating the weights in (|5.2p represents a considerable challenge in 
such a high-dimensional system. Moreover, as already indicated, this expression is derived 
under special conditions not satisfied by the Navier-Stokes equations (e.g. hyperbolicity) 
as well as being derived only for (unstable) periodic orbits rather than (unstable) relative 
periodic orbits (the latter do not exist in very low dimensional systems). A modified 
theory is being developed (Cvitanovic 2012) to include them but here, in the spirit of 
what follows, we brush over this subtlety to treat relative periodic orbits just like periodic 
orbits (see also Lopez et al. 2005). Given these issues, we proceed in a more pragmatic 
fashion in keeping with the previous suggestive work of Zoldi & Greenside (1998) and 
Kazantsev (1998,2001) (see also Dettmann & Morriss 1997). These authors proposed 
and tested a strategy of developing weights based upon only the unstable eigenvalues 



associated with the recurrent flow (i.e. 5fte(A^*^) > 0). In particular, Zoldi 
(1998) advocated the 'escape-time' weighting 

1 



Wi oc 



protocol 1 



Greenside 



(5.3) 



where dC'-*^ is the set of k such that 5fte(A^''') > since this efficiently captures how 
unstable the recurrent flow is and inversely correlates this with how long the turbulent 
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Figure 13: A sampling of the recurrent structures found at Re = 40 plotted over the 
zoom-in box shown in figure El Upper left shows P1,P2,R2 and i?4; Upper right shows 
R28,R29 and i?38; Lower left shows i?19, i?39 and i?55; and lower right shows i?40, i?48 
and R53. In all cases, symbols added to lines are to assist in their distinction and are 
placed 1 time unit apart to indicate speed of flow. 



trajectory should spend in its vicinity. Kazantsev (1998, 2001) argued that this formula 
should be modified to reflect the fact that longer period orbits have a greater 'presence' 
in phase space than shorter period orbits and added the period T*-*^ to the numerator 

Wi oc — protocol 2 (5.4) 

Significantly, this protocol suppresses any contribution from equilibria or travelling waves. 
We now test both of these protocols together with a 'control' choice of 'no weighting' so 
just 

Wi (X 1 protocol 3. (5-5) 
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Figure 14: The scaled, time-averaged dissipation Re{D / Diam)t verses Re for recurrent 
flows discovered by the series A runs. Blue thin hncs trace recurrent flows found at 
Re = 40, red thin dashed hues new recurrent flows found at Re = 60, thick cyan lines 
those found at Re = 80 and very thick magenta lines those at Re ~ 100. The black 
dots identify the subset of solutions which were identified by processing the dns runs at 
the respective Re as opposed to just being continued up or down from other Re (e.g. 
there are 6 dots at Re = 60 corresponding to El, Tl, TS, T4, R7 and R8: R56-R58 were 
discovered in the series B runs, and none of the red dashed lines join dots at Re = 40). 
Open circles indicate limits beyond which a solution branch could not be continued. The 
situation is clearly complicated with solutions seemingly dynamically important at some 
Re but not at others. 

The key measures we use to characterise the 2D turbulence simulated here are pdfs of 
E{t) and D{t) together with the profiles u{y), Urms{y) and Vrms{y) (the pdf of I{t) was 
also considered but adds little information to that provided by the pdf for D{t)). 

5.1. i?e = 40 

The pdfs of E{t) / Eiam and D{t) / Diam at Re = 40 are shown in figure [20] along with 
the predictions using protocols 1-3, that is, the individual pdfs of each recurrent flow 
are weighted together appropriately to produce an overall pdf as in (|5.1[) . In the case of 
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Figure 15: The travelling waves T3 (left), T4 (middle) and T5 (right) at Re = 100. 
Vortieity is contoured using 15 contours between -20 (dark, red) and 15 (light, white) 
(range is -18.9 ^ w 8.03 for T3, -11.5 ^ w < 11.5 for TA and -13.7 ^ w 13.7 for 
T5). 
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Figure 16: The normalised dissipation D{t) / Diam verses I / Diam for the recurrent 
flows found at Re ~ 60 with the pdf of the DNS turbulence plotted in the back- 
ground (11 shades at levels 10" where a = — 5, — 4.5, . . . , — 0.5, 0). Plotted are El 
{D/Diam = I /Diam = 0.091), Tl (0.033), T3 (0.068), T4 (0.034), i?7, i?8, i?56 - 8. 

E{t) I Eiam all three predictions are very good at the central peak of the pdf but each fails 
to capture the clear shoulder at higher energies in the DNS albeit at a level of the pdf 
a factor of ~ 30 lower. There is a similar story for the D{t) / Diam pdf although in this 
case, the performance of the stability-motivated protocols 1 and 2 seem noticeably more 
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Figure 17: The normahsed dissipation D{t)/Diam verses I/Diam for the recurrent 
flows found at Re = 80 with the pdf of the DNS turbulence plotted in the back- 
ground (11 shades at levels 10" where a = — 5, — 4.5, . . . , — 0.5, 0). Plotted are Tl 
{D/Diam = I/Diara = 0.0196), T3 (0.0514), T5 (0.0435, yellow dot), R7 - R13 and 
R59 - R60. 



effective in capturing the DNS pdf. Again, the extremes of the DNS pdf are not captured 
reflecting the fact, as commented earlier, that perhaps not enough recurrent flows with 
large or small dissipation excursions have been found. Given the reduced value of the pdf 
there, these are infrequently visited and thus require very long runs to realistically have 
a chance to extract them. As already mentioned, what initially looked like long runs of 
10^ time units were actually not long enough. As further independent evidence of this, 
plots of the mean profiles u{y) from each 'long' run at the same Re (see Table 1) showed 
noticeable differences between each other and to the expected asymptotic state which 
respects all the symmetries of the system. In particular, the obvious symmetry that the 
mean profile should be invariant under 7r/2 shifts in y was clearly violated. To ameliorate 
this, we decided to 'symmetrise' the DNS mean profile by extracting that part {U^^{y)) 
from the signal (u) which does satisfy all the symmetries listed in §2.11 Explicitly 

[/^''(y) ^ E S-^C^'^CS^y) with U" -.^ -[u{y) + 3i-'uiJiy)] (5.6) 

m— 

(recall n = 4) and similarly for uf,f,;^ and w^ms- This process picks out the following 
Fourier coefficients 

U^'^iy) ■■= E "™ sin(4(2m + l)y) ( v^^, ) = ^ ( ) cos{8my). (5.7) 

m— m— 
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Figure 18: The normalised dissipation D{t)/ Diam verses I/Diam for the recurrent flows 
found at Re = 100 with the pdf of the DNS turbulenee plotted in the background 
(11 shades at levels 10" where a = -5, -4.5, . . . , -0.5, 0). Plotted are Tl (D/Diam = 
I/Diam = 0.0137), T3 (0.0353), TA (0.0170), T5 (0.0297, yellow dot), P4, R7, Rll, R12, 
RU-R17 and R61. 



from the complete Fourier series for u, u„ns and Vms- In particular, the symmetrised 
mean profile U^^{y) has the leading Fourier modal form of sin4j/, which mimicks the 
forcing, and a leading correction of sin 12?/. Such a profile needs only be plotted over 
y G [0, 7r/4] which is done in figure !^ along with the predictions from protocols 1-3. This 
comparison looks impressive with qq — 0.232 (cf expression (|5.7I) ) in the DNS, compared 
to 0.214 (protocol 1), 0.215 (2) and 0.219 (3). For aU, ai = 0(10"^) and 02 = 0(10"^) 
(as way of comparison, the 'raw' mean flow u has a leading non-symmetrised part given 
by 0.0280 cosy + 0.0268 siny, i.e. roughly 10% smaller than the symmetrised part). The 
explanation for why the (symmetrised) mean profile matches the forcing profile so well 
is currently unclear to us and it is tempting to speculate that actually a„ — >■ {n ^ 1) 
with the period of averaging. Sarris et al. (2007) study the statistics of 3D Kolmogorov 
flow for various computational domains and use two measures to signal whether their 
statistics have converged sufficiently over a period of time integration. One, 



{{I)t - {D)t) 



72 := (5-8) 



(equation (22), Sarris et al. 2007) assesses the extent to which the energy input into the 
flow matches the the energy dissipated and is easily calculated from our output data: 72 
is at most O(10~^) for all our 10^ time unit runs. Even with this small value, we find 
evidence that the mean profile is far from converged to what is expected (i.e. satisfies all 
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Figure 19: (Low resolution version for the arXiv) The dns trajectory (left) synchronised 
with the subsequently converged periodic orbit Pi (right) at Re = 100. Time proceeds 
downwards with snapshots at to, to + 0.2, to + 0.6, to + 0.8 and to + 1.0 (the period of P4 
is 1.185). The vorticity scale ranges from -26 (dark red) to 12 (light white). 
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Figure 20: The probability density functions for E{t)/Eiam and D{t)/Diam from DNS 
(blue thick line) and predictions using weighting protocol 1 (red, thick line with dots) , 2 
(green, thin line) and 3 (black thick dashed) at Re = 40. 40 bins were used to calculate 
the pdfs for the recurrent flows and 100 bins for the DNS due to its greater range. These 
choices gave the best balance of resolution with the data available. 



the symmetries of the problem) which emphasizes how easy it is to unwittingly collect 
unconverged statistics. 

Figure [2T] also shows the equivalent plot for M^ms and Wrms- For u„„s, bo = 0.687 for the 
DNS which clearly differs from aU the predictions - 0.272 (protocol 1), 0.289 (2) and 0.276 
(3). The comparison for Vmis, however, is much better: cq — 0.933 for the DNS verses 
0.887 (protocol 1), 0.891 (2) and 0.879 (3). One possible reason why the Urms comparison 
is poor is that Urms is calculated in the DNS 'on the fly' by subtracting the current best 
estimate of the mean u from the current streamwise velocity (see I2.12|) rather than using 
the final mean profile to a posteriori calculate the streamwise fluctuation field. Given our 
realisation now that the mean profile takes a long time to converge to its (symmetric) 
asymptotic state, there is likely to be a significant error (henceforth we consider only 
Vrms for Re = 60, 80 and 100). 

An inescapable conclusion from these comparisons so far is that the 'control' protocol 
3 of actually 'no weighting' performs almost as well as the other stability-motivated 
protocols 1 and 2. It's worthwhile at this point to clarify why. The upper plot in figure 
[22l shows how the peak symmetrised mean value, U^^{tt/8), varies for each recurrent fiow 
and compares these values with the DNS and predictions from the 3 protocols. From this 
it is clear that most of the recurrent flows are good predictors individually and so when 
they are mixed together the result is still reasonably good. The lower plot in figure [H] 
shows how the weights vary across the recurrent flows in the 3 different protocols. Again, 
there is not that much variation over the majority (although notice that Wi ^ for 
i = 1,2,3 in protocol 2 since T^*) = 0) which presumably reflects the fact that the 
stability characteristics of the recurrent flows are all pretty similar. 
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Figure 21: The symmetrised mean flow U^^{y) (left) and symmetrised u^^^{y) (right 
plot left part) and v^^g{y) (right part of right plot) from DNS (blue thick line) and 
predictions using weighting protocol 1 (red, thick line with dots), 2 (green, thin line) and 
3 (black thick dashed) at Re = 40. 

5.2. Re = 60,80 & 100 

The smaller number of recurrent flows extracted for Re > 40 means that it is harder 
to generate reasonably smooth predictions for the pdfs of the energy and dissipation. 
While the same number of 100 bins as at Re = 40 can be used to generate a smooth 
DNS pdf, only 60 bins could be used to sum the pdfs of the recurrent flows. This number 
produced the best compromise of granularity across the range while ensuring that there 
is enough data in each bin for (reasonable) smoothness at least for Re = 60 and Re = 80 
(the sparse coverage of the dissipation range by the recurrent flows at Re = 100 - see 
figure [TH] - prevented any useful plot to be generated) . Figure [23] shows the result of this 
procedure for the dissipation pdf at Re = 60 and 80 (plots not shown for the energy 
are similar). Here, protocol 2 offers the best partial fit both at Re = 60 and Re = 80 
although it is clear that much of the turbulent dissipation range extends higher than 
any of the recurrent fiows found at both Re (also clear from \W\ and [T7|) and that the 
'prediction' is of limited quality. It's also worth remarking that the predictions are now 
more distinquishable at these higher Re which can be traced back to the separating 
stability characteristics of the recurrent flows. For example. El has J^^^i^j) — 5.053 
whereas this is just 0.139 for Tl. 

Figure [24] shows the symmetrised mean profile as calculated from the DNS and pre- 
dicted by the 3 protocols at Re = 60 and Re = 80. Again, somewhat paradoxically, the 
'control' protocol does the best job in both cases. At Re = 60 oq = 0.2277, 0.2298 and 
0.2296 (cf expression ()5.7|) ) across the three series A DNS runs listed in Table 1 (with 
again « 10% non-symmetrised part in all cases). In comparison, gq = 0.148 (protocol 
1), 0.175 (2) and 0.200 (3). For Vrms, cq = 1.099 in the DNS run to be compared with 
the predictions of 1.064 (protocol 1), 1.057 (2) and 1.055 (3). At Re = 80, oq = 0.2009, 
0.2016 and 0.1998 across the three series A DNS runs listed in Table 1. In comparison, 
ao = 0.1416 (protocol 1), 0.1403 (2) and 0.1944 (3). 
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Figure 22: Top: The peak symmetrised mean U^^{tt/8) for the DNS (thick blue horizontal 
line), prediction using protocol 1 (red dotted line), prediction using protocol 2 (green 
line), the control prediction (dashed black line) and for each recurrent flow listed in the 
order given in Table 1 at Re = 40 (so i = 20 is i?25 and i = 45 is i?50 marked with black 
dotted lines). Bottom: the weights for the three protocols (1-red line with dots, 2-green 
line, 3-black horizontal line) plotted for all the recurrent structures at Re ~ 40. Again 
the vertical (black) dotted lines indicate i?25 and R50. 



Finally, we note that at Re = 100, oq = 0.1777, 0.1782 and 0.1783 across the three 
series A DNS runs listed in Table 1 (with now w 20% non-symmetrised part in all cases) . 
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Figure 23: The probability density functions for D[t)/Diam from DNS (blue thick line, 
cyan and magenta lines) and predictions using weighting protocol 1 (red, thick line with 
dots), 2 (green, thin line with triangles), and 3 (black thick dashed with squares) at 
Re = 60 (left) and Re = 80 (right). 60 bins were used to calculate the pdfs for the 
recurrent flows and 100 bins for the DNS due to its greater range. These choices gave 
the best balance of resolution with the data available. 
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Figure 24: The symmetrised mean flow U^^'{y) from DNS (blue thick line) and predictions 
using weighting protocol 1 (red, thick line with dots), 2 (green, thin line) and 3 (black 
thick dashed) for Re = 60 (left) and Re ^ 80 (right). 
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6. Discussion 

In this study, we have considered 2D turbulence on the torus [0, 27r]^ forced nionochro- 
matically in one direction (Kolmogorov flow) . By looking for near recurrences of the flow 
in long direct numerical simulation runs, sets of exactly closed flow solutions 'embedded' 
in this turbulence have been extracted at different forcing amplitudes (Re). We have then 
tried to use these sets of recurrent flows to reconstruct key statistics of the turbulence 
motivated by Periodic Orbit Theory in low-dimensional chaos. The approach has been 
reasonably successful at Re = 40 - see figures [201 and [2T] - where 50 recurrent flows were 
found with the majority buried in the part of phase space most populated by the tur- 
bulence. In contrast, at Re ~ 60, 80 and 100, the limited size of the recurrent flow sets 
found has made the approach largely impotent. Even at Re = 40, the success achieved 
seems more reliant on just extracting lots of similar-looking recurrent flows buried in the 
most popular part of phase space for the turbulence than on any sophisticated choice 
of weighting coefficients. Indeed, one is reminded of Kawahara & Kida's (2001) conclu- 
sion that one judiciously chosen periodic orbit is 'enough' to be a valuable proxy of the 
turbulence. We can sympathise with this viewpoint but only if the comparison with the 
turbulence statistics is not too demanding. The key issue, of course, plaguing this in- 
vestigation is the paucity of recurrent flows found from the finite DNS data generated. 
This is perhaps the main message to come out of this work: Periodic Orbit theory for 
fluid turbulence is a promising approach but only if enough - O(IOO)? - recurrent flows 
are gathered which requires very long turbulence data sequences. A time sequence of 10^ 
time units seems marginally adequate for Re = 40 but is maybe two orders of magnitude 
too short for Re = 60 and beyond. Unfortunately, without these large sets, it has been 
impossible to discern between different weighting protocols. 

Operationally, the work described here has been time-consuming both computation- 
ally in generating near-recurrence episodes and attempting to converge them, as well 
as 'manually' because of all the careful processing (e.g. calculating their stability) and 
archiving of the recurrent flows needed (e.g. does a new convergence from a DNS guess 
represent a new recurrent flow or a repeat of a previously extracted flow?). Fortunately, 
there is no reason why this process could not be automated with the objective being 
to 'automatically' generate a basis set of recurrent flows for each Re. Indeed, one could 
hope that such a set at given Re = Rei could be used to predict the turbulent statistics 
at another Re = Re2- This would require each recurrent flow at Rei being continued 
to Re2 and the fresh weights for an expansion being generated from the (new) stability 
information for each recurrent flow - again painstaking work but readily automated. One 
fly in the ointment is the possibility of bifurcations in the interval [Rei, Re2\, particularly 
saddle node bifurcations where two recurrent flows at Rei merge and annihilate before 
Re reaches Re2- Working with large enough recurrent flow sets would presumably smooth 
over this effect somewhat but will not eliminate it entirely. 

Leaving aside these issues for a moment, it's worth re-emphasizing that any recurrent 
flow extracted from DNS data is a simple invariant solution 'buried' in the turbulence. As 
such, each represents a sustained sequence of dynamical processes which contributes to, 
if not underpins, the turbulence itself. Since they are closed in time, they can be analysed 
relatively easily in whatever detail is required to understand key dynamical relationships 
in the flow. This seems a very promising byproduct of the analysis whether one believes 
a Periodic Theory-type expansion of turbulence is possible or not (pursuing this has not 
been the focus here due to the 2-dimensionality of the flow) . 

Finally, the ever-improving computational resources available now have only recently 
made this type of study possible. Even with these, we have underestimated the demands 
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of data collection in 2D turbulence over the small torus [0, 2tt]'^. Major challenges ahead 
include treating large aspect ratio domains - can we find localised recurrent flows? - and 
handling fully 3 dimensional flows - with automated machinery, will the approach be 
practical? There is plenty to explore. 
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Appendix A 

The Ncwton-GMRES-Hook-step algorithm described in the main text is easily ex- 
tended to continue solutions over parameter space such as Re or the domain geometry 
(e.g. a). We briefly describe this extension for solution branch continuation in Re which 
was used to generate figure [T31 A simple strategy is to use the solution x(i?e) as an 
initial guess in the Newton-GMRES-Hook-step algorithm with the hope of converging to 
x(i?e -f (Si?e). This should work provided that 5Re is 'small enough' but is ill-equipped to 
negotiate turning points in the solution branch. A standard, more sophisticated approach 
is arc-length continuation which uses the branch arc-length as a natural, monotonically- 
incrcasing, parametrisation of the solution branch. The key idea is to take small con- 
trollable steps in the arc-length rather than Re. As a result the state vector needs to be 
extended as follows 

X = 

and an extra equation 

dr 



n 

s 
T 
Re 



9x ^ 

dr 



(6.1) 



(6.2) 
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added to determine Re. Previous converged solutions a;(r_i) and x{ro) indicate a rea- 
sonable step size in r, 5r = ro — r_i = ^ {x{rQ) — a;(r_i))2 and allow a prediction to be 
made for the next solution 



dx 

x{ri) « x{ra) + 5r— 
or 



(6.3) 



Given x(ro) and dr, the extra constraint for the Newton method comes from approxi- 
mating (|6.2p as follows 



J{{x'') :-- 



dx 
dr 



(a;" - x{ro)) - (5r w 



for the nth iterate to estimate x{ri). Writing (5a;" := x'^'^^ — a;", then setting 

N(a;"+i) = (5.t".^ |,„ + K(a;") = 
or 

puts the required extra constraint on the iterative improvement 5x"^ . The Newton prob- 
lem p.l3p then becomes 



(6.4) 



(6.5) 













dT 


dRe 


••• (T.flo)^ ••• 









dt 


















F(f2o,so,7o;m) 




5s 









ST 









. SRe _ 




l^ixo) 



dxp ^ 
dr 

(6.6) 

Depending on how easily convergence is obtained, dr can be increased or decreased if 
the algorithm shows signs of divergence. A second order approach to estimating dx/dr 
was actually adopted for the predictive step but the first order estimate proved sufficient 
for the constraint present in (|6.6p . 



